Starting to look at linear regression for predicting pitch pine presence/abundance. I think markdown scripts are easier to keep neat and organize.

All counts are included in this script (PIRI, Shrub Oak, and Other categories)

Libraries:

# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Functions:

# functions -------------------

Import small seedling data:

# global variables -------------------
source("Scripts/SS_Import.R")

Pivot wider to create dataframe where each row is for one plot and has total details for each species group

ss_merge2 <- ss_merge %>% 
  select(Plot_No, Region, Treat_Type, Site, Species_Groups, Total) #this is dropping stump sprout, browse, and germinate data

ss_merge2 <- ss_merge2 %>% 
  pivot_wider(names_from = Species_Groups, values_from = Total)

Import time since data and add it to the small seedling dataset

time_since <- read_csv("CleanData/Treat_Year_Data.csv")

ss_merge3 <- merge(ss_merge2, time_since, by = "Site")
#log transform time from treatment data
ss_merge3$l.TFT <- log(ss_merge3$Time_from_Treat)

Run the ‘Add_BA’ script and merge with dataset:

source("Scripts/Add_BA.R")

# merge with ss dataset -------------------
ss_merge4 <- merge(ss_merge3, prism_BA, by = "Plot_No")

Run ‘Ground_Data.R’ script and add it to small seedling dataset:

source("Scripts/Ground_Data.R")

# merge with ss dataset -------------------
ss_merge5 <- merge(ss_merge4, ground3, by = "Plot_No")

rm(ss_merge2,
   ss_merge3,
   ss_merge4,
   time_since)

Source and add veg data

source("Scripts/Veg_Data.R")

# merge with ss dataset 
ss.all <- merge(ss_merge5, veg3, by = "Plot_No")

Create log transformed categories for newly added variables, then select for just the desired variables:

ss.all$l.PIRI <- log(ss.all$PIRI + 1)
ss.all$l.SO <- log(ss.all$Shrub_Oak + 1)
ss.all$l.other <- log(ss.all$Other + 1)

ss.all2 <- ss.all %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)

Select just for numerical vs log and then look at paired plots:

#not transformed
ss.num <- ss.all2 %>% 
  select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(ss.num, aes(color = Treat_Type))
ggpairs(ss.num)

#log transformed
ss.numl <- ss.all2 %>% 
  select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)

ggpairs(ss.numl, aes(color = Treat_Type))
ggpairs(ss.numl)

rm(ss.num,
   ss.numl)

Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)

The above script and naming conventions are the same as the SS_GLMM script

Full dataset is called ss.all2.

Hidden below is the original analysis and step by step model elimination conducted without treatment type included. Best model without treatment type had a higher AIC and is underdispersed according to DHARMa tests

Pitch pine GLMM

Model step wise elimination with treatment type included

ss.m1 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m1) #724.86
## [1] 724.8561
#would like to test ba vs ba piri
ss.m2a <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m2a) #723.6
## [1] 723.6111
lrtest(ss.m1, ss.m2a) # p = .4
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + 
##     l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  14 -348.43                     
## 2  13 -348.81 -1 0.7551     0.3849
ss.m2b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m2b) #725.9
## [1] 725.919
lrtest(ss.m1, ss.m2b) # p =.08, so not important
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  14 -348.43                       
## 2  13 -349.96 -1 3.0629     0.0801 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m2a, ss.m2b)

# test ld vs mineral exposure
ss.m3a <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m3a) #723.0
## [1] 722.9976
lrtest(ss.m1, ss.m3a) #p=.25
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  14 -348.43                     
## 2  11 -350.50 -3 4.1415     0.2466
ss.m3b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m3b) #730.1 
## [1] 730.1247
lrtest(ss.m1, ss.m3b) # p = 0.01, loss of avg ld significant
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  14 -348.43                       
## 2  11 -354.06 -3 11.269    0.01036 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m1, ss.m3b)

# test l.other
ss.m4 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m4) #722.6
## [1] 722.5895
lrtest(ss.m3a, ss.m4) #p=.2, can drop other
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  11 -350.50                    
## 2  10 -351.29 -1 1.592      0.207
rm(ss.m3a)

# test l.SO
ss.m5 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m5) #725.5
## [1] 725.4735
lrtest(ss.m4, ss.m5) #p = 0.03, l.so important
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | 
##     Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)  
## 1  10 -351.29                      
## 2   9 -353.74 -1 4.884    0.02711 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m5)

#test veg
ss.m6 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(ss.m6) #727
## [1] 726.9952
lrtest(ss.m4, ss.m6) # p = 0.011, veg important
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  10 -351.29                       
## 2   9 -354.50 -1 6.4056    0.01138 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m6)

do pairwise on variables of interest from reduced model

ss.pair <- ss.all2 %>% 
  select(Region, Treat_Type, l.PIRI, l.SO, l.Veg_Total, avgLD_l, l.TFT)

ggpairs(ss.pair, aes(color = Treat_Type))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

rm(ss.pair)

Not seeing any clear, strong correlations

Best model:

ss.m4 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
summary(ss.m4) #722.6
##  Family: poisson  ( log )
## Formula:          
## PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +  
##     (1 | Site/Plot_No)
## Data: ss.all2
## 
##      AIC      BIC   logLik deviance df.resid 
##    722.6    764.7   -351.3    702.6      488 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance Std.Dev.
##  Plot_No:Site (Intercept) 4.657    2.158   
##  Site         (Intercept) 2.247    1.499   
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Conditional model:
##                    Estimate Std. Error z value   Pr(>|z|)    
## (Intercept)         -5.2307     1.7872  -2.927    0.00343 ** 
## Treat_TypeFallRx     6.8811     1.4362   4.791 0.00000166 ***
## Treat_TypeHarvest    6.4777     1.5012   4.315 0.00001595 ***
## Treat_TypeMowRx      5.2753     1.4183   3.719    0.00020 ***
## Treat_TypeSpringRx   6.7837     1.4407   4.709 0.00000249 ***
## l.SO                -0.5172     0.2390  -2.164    0.03049 *  
## avgLD_l             -1.1967     0.4340  -2.757    0.00583 ** 
## l.Veg_Total         -0.8041     0.3280  -2.451    0.01423 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#This is exactly the same, just checking Maria's comment
ss.m4a <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site) + (1|Plot_No),
                 data = ss.all2,
                 family = poisson)
summary(ss.m4a)
##  Family: poisson  ( log )
## Formula:          
## PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +  
##     (1 | Site) + (1 | Plot_No)
## Data: ss.all2
## 
##      AIC      BIC   logLik deviance df.resid 
##    722.6    764.7   -351.3    702.6      488 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name        Variance Std.Dev.
##  Site    (Intercept) 2.247    1.499   
##  Plot_No (Intercept) 4.657    2.158   
## Number of obs: 498, groups:  Site, 47; Plot_No, 498
## 
## Conditional model:
##                    Estimate Std. Error z value   Pr(>|z|)    
## (Intercept)         -5.2307     1.7872  -2.927    0.00343 ** 
## Treat_TypeFallRx     6.8811     1.4362   4.791 0.00000166 ***
## Treat_TypeHarvest    6.4777     1.5012   4.315 0.00001595 ***
## Treat_TypeMowRx      5.2753     1.4183   3.719    0.00020 ***
## Treat_TypeSpringRx   6.7837     1.4407   4.709 0.00000249 ***
## l.SO                -0.5172     0.2390  -2.164    0.03049 *  
## avgLD_l             -1.1967     0.4340  -2.757    0.00583 ** 
## l.Veg_Total         -0.8041     0.3280  -2.451    0.01423 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Checking model fit:

ss.m4_sr <- simulateResiduals(ss.m4, n = 1000, plot = TRUE) #looks pretty good

testResiduals(ss.m4_sr) #passes uniformity, dispersion, and outliers

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.055502, p-value = 0.09301
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0070259, p-value = 0.212
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00118473895582329 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.055502, p-value = 0.09301
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0070259, p-value = 0.212
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00118473895582329 ) 
##                                                  0
testQuantiles(ss.m4_sr) #passes

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.2776
## alternative hypothesis: both
testZeroInflation(ss.m4_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99615, p-value = 0.946
## alternative hypothesis: two.sided

Model 4 has the lowest AIC and great model fit

Time to pairwise test

emmeans(ss.m4, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type     rate       SE  df asymp.LCL asymp.UCL
##  Control    0.000147 0.000196 Inf 0.0000108     0.002
##  FallRx     0.143092 0.098074 Inf 0.0373430     0.548
##  Harvest    0.095594 0.081151 Inf 0.0181066     0.505
##  MowRx      0.028723 0.019958 Inf 0.0073586     0.112
##  SpringRx   0.129822 0.093714 Inf 0.0315422     0.534
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast             ratio      SE  df null z.ratio p.value
##  Control / FallRx   0.00103 0.00147 Inf    1  -4.791  <.0001
##  Control / Harvest  0.00154 0.00231 Inf    1  -4.315  0.0002
##  Control / MowRx    0.00512 0.00726 Inf    1  -3.719  0.0019
##  Control / SpringRx 0.00113 0.00163 Inf    1  -4.709  <.0001
##  FallRx / Harvest   1.49687 1.48518 Inf    1   0.407  0.9943
##  FallRx / MowRx     4.98172 4.23482 Inf    1   1.889  0.3232
##  FallRx / SpringRx  1.10222 0.97175 Inf    1   0.110  1.0000
##  Harvest / MowRx    3.32809 3.28617 Inf    1   1.218  0.7411
##  Harvest / SpringRx 0.73635 0.73882 Inf    1  -0.305  0.9981
##  MowRx / SpringRx   0.22125 0.18896 Inf    1  -1.766  0.3935
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

FallRx, Harvest, MowRx, SpringRx all significantly different from Control - no treatment significantly different from another

Model ss.m4 uses Treatment Type, Scrub Oak presence, average leaf litter depth, and total veg cover as predictors. Make a graph for each

I’d like to see some graphs before I close out

library(sjPlot)
library(sjlabelled)
## 
## Attaching package: 'sjlabelled'
## The following object is masked from 'package:forcats':
## 
##     as_factor
## The following object is masked from 'package:ggplot2':
## 
##     as_label
## The following object is masked from 'package:dplyr':
## 
##     as_label
library(sjmisc)
## 
## Attaching package: 'sjmisc'
## The following object is masked from 'package:purrr':
## 
##     is_empty
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## The following object is masked from 'package:tibble':
## 
##     add_case
set_theme(base = theme_classic(),
          theme.font = 'serif',
          axis.title.size = 1.5,
          axis.textsize.x = 1.5,
          axis.textsize.y = 1.5,
          title.size = 2.5,
          title.align = "center",
          legend.pos = "right",
          legend.size = 1.5,
          legend.title.size = 1.5,
          #legend.bordercol = "black",
          legend.item.size = .75)

plot_model(ss.m4) #this is incidence rate ratios

plot_model(ss.m4, type = "diag") #random vs normal quantiles
## $`Plot_No:Site`
## `geom_smooth()` using formula = 'y ~ x'

## 
## $Site
## `geom_smooth()` using formula = 'y ~ x'

plot_model(ss.m4, type = "re") #this plots random effects
## [[1]]

## 
## [[2]]

plot_model(ss.m4, 
           type = 'pred', 
           terms = 'Treat_Type') #plot marginal effects (i might have done this wrong)
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

# SS PIRI plot 1 ------------------- LD vs. TT
plot_model(ss.m4, 
                       type = 'pred', 
                       terms = c('avgLD_l', 'Treat_Type'),
                       axis.title = c("Average Leaf Litter Depth (log transformed)", "Total Count of Pitch Pine"),
                       title = "Predicted Counts of Pitch Pine Seedlings <50cm",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

# SS PIRI plot 2 ------------------- Veg vs TT
plot_model(ss.m4, 
                       type = 'pred', 
                       terms = c('l.Veg_Total', 'Treat_Type'),
                       axis.title = c("Average Understory Vegetation Cover (log transformed)", "Total Count of Pitch Pine"),
                       title = "Predicted Counts of Pitch Pine Seedlings <50cm",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

# SS PIRI pot 3 ------------------- Shrub Oak seedling vs TT
plot_model(ss.m4, 
                       type = 'pred', 
                       terms = c('l.SO', 'Treat_Type'),
                       axis.title = c("Shrub oak small seedling counts (log transformed)", "Total Count of Pitch Pine"),
                       title = "Predicted Counts of Pitch Pine Seedlings <50cm",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

# SS PIRI plot 4 ------------------- Time from Treatment vs TT
plot_model(ss.m4, 
                       type = 'pred', 
                       terms = c('l.TFT', 'Treat_Type'),
                       axis.title = c("Time from Treatment (log transformed)", "Total Count of Pitch Pine"),
                       title = "Predicted Counts of Pitch Pine Seedlings <50cm",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

Next steps:

2e. other interaction terms to think about? 3. take another look at treat type glmms 4. looks at lr sans treat type. Better AIC or fit? stronger predictors? maybe in LS or SA models? 6. start with full models post scatter plot analysis 7. are there other ways that I want to think about veg cover? e.g., just pulling out SO cover? VA/GA cover?

Shrub oak GLMM

Begin model elimination

so.ss1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss1) #1897.4
## [1] 1897.374
#test piri BA
so.ss2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss2) #1896
## [1] 1895.956
# test total ba
so.ss3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss3) #1894.1
## [1] 1894.07
# test mineral
so.ss4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss4) #1892.3
## [1] 1892.281
# test litter depth
so.ss5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss5) #1890.1
## [1] 1890.831
# test other seedling
so.ss6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss6) #1889.5
## [1] 1889.516
lrtest(so.ss5, so.ss6) # p = .4
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  10 -935.42                    
## 2   9 -935.76 -1 0.685     0.4079
# test PIRI seedling
so.ss7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss7) #1892.6
## [1] 1892.648
lrtest(so.ss6, so.ss7) # p = 0.02, so l.PIRI is important
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1   9 -935.76                       
## 2   8 -938.32 -1 5.1315     0.0235 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# add PIRI back in, test veg cover
so.ss8 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss8) #1888.9
## [1] 1888.943
lrtest(so.ss6, so.ss8) # p = 0.2
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   9 -935.76                     
## 2   8 -936.47 -1 1.4263     0.2324
# test treat type
so.ss9 <- glmmTMB(Shrub_Oak ~ l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
AIC(so.ss9) #1911.5
## [1] 1911.528
lrtest(so.ss9, so.ss8) # p > 0.001, so TT is important
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq  Pr(>Chisq)    
## 1   4 -951.76                          
## 2   8 -936.47  4 30.585 0.000003719 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Best model, then test

so.ss8 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = poisson)
summary(so.ss8)
##  Family: poisson  ( log )
## Formula:          
## Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ss.all2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1888.9   1922.6   -936.5   1872.9      490 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance Std.Dev.
##  Plot_No:Site (Intercept) 1.923    1.387   
##  Site         (Intercept) 1.485    1.219   
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Conditional model:
##                    Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)         -4.3474     0.4346 -10.003 < 0.0000000000000002 ***
## Treat_TypeFallRx     3.1111     0.6183   5.032         0.0000004851 ***
## Treat_TypeHarvest    1.3851     0.7377   1.877               0.0605 .  
## Treat_TypeMowRx      3.1769     0.5759   5.516         0.0000000346 ***
## Treat_TypeSpringRx   2.9035     0.6316   4.597         0.0000042769 ***
## l.PIRI              -0.4518     0.1890  -2.390               0.0168 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ss8_sr <- simulateResiduals(so.ss8, n = 1000, plot = TRUE) #passes KS, quantile deviations signficant

testResiduals(so.ss8_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.033251, p-value = 0.6406
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0083534, p-value = 0.002
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.96
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00184738955823293 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.033251, p-value = 0.6406
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0083534, p-value = 0.002
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.96
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00184738955823293 ) 
##                                                  0
testDispersion(so.ss8_sr, alternative = "less") #fails, so it is underdispersed

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0083534, p-value = 0.001
## alternative hypothesis: less
#try again with gen pois / com pois ditros


so.ss10 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = compois)
summary(so.ss10)
##  Family: compois  ( log )
## Formula:          
## Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ss.all2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1869.6   1907.5   -925.8   1851.6      489 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance Std.Dev.
##  Plot_No:Site (Intercept) 1.191    1.091   
##  Site         (Intercept) 1.504    1.226   
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Dispersion parameter for compois family (): 4.39 
## 
## Conditional model:
##                    Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)         -4.0292     0.4373  -9.213 < 0.0000000000000002 ***
## Treat_TypeFallRx     3.0398     0.6152   4.941         0.0000007773 ***
## Treat_TypeHarvest    1.4196     0.7350   1.931               0.0534 .  
## Treat_TypeMowRx      3.1285     0.5737   5.453         0.0000000494 ***
## Treat_TypeSpringRx   2.9044     0.6288   4.619         0.0000038500 ***
## l.PIRI              -0.4238     0.1779  -2.382               0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ss10_sr <- simulateResiduals(so.ss10, n = 1000, plot = TRUE)

testResiduals(so.ss10_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.040131, p-value = 0.3989
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.050676, p-value = 0.01
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00146586345381526 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.040131, p-value = 0.3989
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.050676, p-value = 0.01
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00146586345381526 ) 
##                                                  0
testDispersion(so.ss10_sr, alternative = "less") # still under dispersed

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.050676, p-value = 0.005
## alternative hypothesis: less
so.ss11 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
                 data = ss.all2,
                 family = genpois)
summary(so.ss11)
##  Family: genpois  ( log )
## Formula:          
## Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ss.all2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1854.2   1892.1   -918.1   1836.2      489 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance Std.Dev.
##  Plot_No:Site (Intercept) 0.2882   0.5369  
##  Site         (Intercept) 0.9904   0.9952  
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Dispersion parameter for genpois family (): 5.87 
## 
## Conditional model:
##                    Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)         -3.4252     0.3608  -9.493 < 0.0000000000000002 ***
## Treat_TypeFallRx     2.8087     0.5032   5.581       0.000000023864 ***
## Treat_TypeHarvest    1.4256     0.6233   2.287               0.0222 *  
## Treat_TypeMowRx      2.9454     0.4699   6.268       0.000000000366 ***
## Treat_TypeSpringRx   2.8173     0.5213   5.404       0.000000065125 ***
## l.PIRI              -0.3657     0.1623  -2.253               0.0242 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ss11_sr <- simulateResiduals(so.ss11, n = 1000, plot = TRUE)

testResiduals(so.ss11_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.025101, p-value = 0.9123
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.28065, p-value = 0.092
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00166666666666667 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.025101, p-value = 0.9123
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.28065, p-value = 0.092
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00166666666666667 ) 
##                                                  0
testDispersion(so.ss11_sr, alternative = "less") # at 0.05, so could possibly work

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.28065, p-value = 0.046
## alternative hypothesis: less
testQuantiles(so.ss11_sr) #quantiles are not good

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.00001043
## alternative hypothesis: both
emmeans(so.ss11, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type response     SE  df asymp.LCL asymp.UCL
##  Control       0.183 0.0666 Inf    0.0898     0.373
##  FallRx        3.038 1.1048 Inf    1.4896     6.196
##  Harvest       0.762 0.4032 Inf    0.2701     2.149
##  MowRx         3.483 1.1174 Inf    1.8574     6.532
##  SpringRx      3.064 1.2098 Inf    1.4134     6.643
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast            ratio     SE  df null z.ratio p.value
##  Control / FallRx   0.0603 0.0303 Inf    1  -5.581  <.0001
##  Control / Harvest  0.2404 0.1498 Inf    1  -2.287  0.1489
##  Control / MowRx    0.0526 0.0247 Inf    1  -6.268  <.0001
##  Control / SpringRx 0.0598 0.0312 Inf    1  -5.404  <.0001
##  FallRx / Harvest   3.9873 2.5162 Inf    1   2.192  0.1827
##  FallRx / MowRx     0.8722 0.4141 Inf    1  -0.288  0.9985
##  FallRx / SpringRx  0.9914 0.5212 Inf    1  -0.016  1.0000
##  Harvest / MowRx    0.2188 0.1325 Inf    1  -2.509  0.0885
##  Harvest / SpringRx 0.2486 0.1605 Inf    1  -2.156  0.1966
##  MowRx / SpringRx   1.1367 0.5632 Inf    1   0.259  0.9990
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

FallRx, MowRx, and SpringRx significantly different than control; none of the treatments significantly different from each other

Doing some data exploration via scatterplots

This analysis can be seen in the pair plots at the top of script

Veg cover, average litter depth, and time from treatment all have a clear negative relationship (as they decrease, total increases), some more significant than others

Basal area per hectare and mineral soil exposure have positive relationship (as they increase, total increase), these relationships appear weaker than the above negative relationships

More scatterplots hidden below; see pairwise at beginning of script

Hidden below are previous, discarded PIRI models :